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Introduction 

This  project  was  concerned  with  the  development  of  mathematical  methods  and  com¬ 
putational  procedures  for  reliable  quantitative  determination  of  margins  of  safety  (MS)  for 
aircraft  structural  systems.  By  definition: 

MS  =  -  1  (1) 

0.MAX 

where  <2>all  >  0  is  the  allowable  value  of  a  functional  prescribed  by  design  specifications 
and  </>max  >  0  is  the  value  of  that  functional  for  a  given  loading  condition.  The  principal 
objective  of  this  project  was  to  develop  procedures  for  the  estimation  and  control  the  errors 
associated  with  numerical  determination  of  <A\ iax  and  proper  interpretation  of  experiments 
aimed  at  the  definition  of  0all-  The  investigation  focused  on  two  areas:  Failure  through 
loss  of  stability  in  composite  plates  and  shells  and  failure  of  adhesively  bonded  joints  in 
composite  materials  subjected  to  thermal  and  mechanical  loads.  There  are  many  important 
practical  applications,  such  as  electronic  packaging,  ceramics,  adhesively  bonded  joints  and 
laminated  composites. 

The  specific  objectives  of  the  project  were: 

1.  Investigate  p- adaptive  procedures.  The  goal  was  to  identify  specific  regions  on  the 
solution  domain  where  a  trial  discretization  needs  to  be  modified. 

2.  Develop  model-adaptive  procedures  for  structural  plates  and  shells  made  of  laminated 
composites. 


20010426  054 


page  2/19 


3.  Develop  methods  and  procedures  for  the  numerical  determination  of  the  eigenpairs 
that  characterize  the  generalized  stress  intensity  factors  along  re-entrant  edges  in  het¬ 
erogeneous  solid  bodies. 

4.  Develop  a  method  for  numerical  determination  of  the  generalized  stress  intensity  factors 
in  heterogeneous  bodies  subjected  to  thermal  and  mechanical  loading. 

5.  Develop  procedures  for  the  simulation  of  time-dependent  material  response  (creep  and 
relaxation). 

One  of  the  difficulties  in  creating  reliable  design  criteria  for  structures  made  of  composite 
materials  is  that  composite  materials  can  fail  by  a  variety  of  mechanisms.  At  present  it  is 
not  possible  to  predict  with  a  reasonable  degree  of  certainty  and  consistency  when  and  how 
a  structure  made  of  composite  materials  will  fail.  Failure  may  be  caused  by  delamination, 
separation  of  the  plies,  buckling  of  the  fibers,  failure  in  the  adhesive,  loss  of  structural 
stability,  or  some  combination  of  these  mechanisms.  Owing  to  the  complexity  of  the  problem, 
estimation  of  margins  of  safety  by  numerical  means  is  not  within  the  state  of  the  art  at 
present.  Progress  is  possible  through  carefully  designed  physical  experiments  coupled  with 
numerical  solution  of  the  models  that  represent  the  experiments.  The  numerical  solutions 
must  have  a  guaranteed  accuracy.  The  development  of  reliable  methods  of  analysis  was  a 
central  consideration  in  this  work.  The  main  points  are  described  in  the  following  section. 

Allowables 

The  definition  and  determination  of  allowable  values  </>all  in  eq.  (1)  requires  close  correlation 
between  experimental  and  computed  data.  The  procedure  involves  the  formulation  of  a 
hypothesis  that  failure  occurs  when  a  certain  functional  <p  reaches  a  critical  value.  The 
functional  cannot  be  observed  directly,- therefore  it  must  be  inferred  from  correlations 
between  computed  values  and  experimental  observations. 

Let  Yij  be  the  zth  ideal  observation  of  the  jth  experiment  and  let  <pi(uEX)  be  the  corre¬ 
sponding  functional  computed  from  the  exact  solution  vfpX  so  that  if  there  were  no  experi¬ 
mental  errors  and  the  mathematical  model  and  the  hypothesis  were  correct  then  we  would 
have 

Yij  -  M*fx)  = 

Due  to  experimental  errors  etj  and  inherent  uncertainties  we  actually  observe  yrj  =  YtJ  ±  et] 
and  compare  this  with  (pi{upE)  where  UpE  is  the  finite  element  solution  corresponding  to  the 
jth  experiment.  Writing: 

Yij  —  =  ^  eij  "Feq  —  0i(u£x)  +  ^i^Fe)  —  FE ) 

v  - - * 

yij 
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and  using  the  triangle  inequality,  we  have 


true  error  apparent  error  error  of  discretization 


This  result  shows  that  to  test  a  particular  hypothesis,  it  is  essential  to  have  both  the  exper¬ 
imental  errors  |e;j|  the  discretization  errors  \4>i(u  FE )  ~  &(uEx)\  under  control,  otherwise  it 
will  not  be  possible  to  know  whether  the  apparent  error  is  dup  to  an  error  in  the  hypothesis 
being  tested;  errors  in  the  numerical  approximation,  or  errors  in  the  experiment.  Further¬ 
more,  control  of  discretization  errors,  in  terms  of  the  data  of  interest,  must  be  intrinsic 
to  the  numerical  solution  method.  The  discretization  errors  must  be  much  less  than  the 
experimental  errors. 

Inherent  uncertainties  in  the  data  may  be  a  dominant  factor.  In  that  case  the  aim  of  the 
experiments  should  include  the  development  of  reliable  statistical  information. 

Summary  of  accomplishments 

The  results  of  this  research  have  been  published,  or  are  being  processed  for  publication. 
A  summary  of  the  work  completed  under  this  project  is  presented  in  the  following. 

1.  p- Adaptive  methods.  In  the  application  of  the  p- version  of  the  finite  element  method 
the  finite  element  mesh  is  fixed  and  the  polynomial  degree  of  the  elements  is  increased 
until  a  desired  level  of  precision  is  reached.  Ideally,  the  mesh  should  be  designed  taking 
into  account  a  priori  information  concerning  the  smoothness  of  the  solution.  Uniform 
meshes  are  optimal  where  the  solution  is  smooth.  In  the  neighborhood  of  singularities 
geometrically  graded  elements  are  optimal  and  in  dimensionally  reduced  problems, 
such  as  plates  and  shells,  the  mesh  design  should  account  for  boundary  layer  effects  as 
well.  In  the  case  of  large  problems  it  is  not  feasible  in  general  to  create  meshes  that 
account  for  singularities  everywhere.  The  goal  is  to  identify  critical  regions,  called 
the  regions  of  primary  interest,  where  stresses  or  stress  intensity  factors  have  to  be 
computed.  In  many  cases  the  regions  of  primary  interest  are  known  a  priori.  In  the 
other  parts  of  the  domain  only  a  reasonably  accurate  representation  of  the  compliances 
is  required.  These  parts  are  called  the  regions  of  secondary  interest.  Compliances  are 
well  approximated  when  the  error  in  strain  energy  is  small. 

The  goal  of  adaptivity  is  to  ensure  that  the  error  in  energy  norm  is  small  and  the  error 
is  nearly  equidistributed  among  the  finite  elements. 

Many  finite  element  analysis  codes,  including  the  three  largest  commercial  codes,  now 
offer  p-extension  capabilities.  In  most  industrial  applications  p-extension  is  used  in 
conjunction  with  mesh  generators  that  typically  generate  nearly  uniform  meshes  of 
tetrahedral  elements.  A  capability  for  adaptive  p-distribution  is  needed  to  ensure 
nearly  optimal  performance  for  a  given  number  of  degrees  of  freedom. 
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Existing  a  posteriori  error  estimators  can  be  grouped  into  three  large  categories:  esti¬ 
mators  based  on  (a)  residuals,  (b)  stress  smoothing,  and  (c)  extrapolation.  Estimators 
of  type  (a)  and  (b)  are  used  in  conjunction  with  the  h-  and  hp- versions,  whereas  esti¬ 
mators  of  type  (c)  require  information  generated  by  p-extensions. 

Construction  of  local  error  indicators  based  on  element-residuals  as  well  as  stress 
smoothing  and  patch  recovery  have  been  the  focus  of  intensive  research.  Indicators 
based  on  residuals  typically  require  local  equilibration  of  the  residuals  and  jumps  in 
tractions  at  inter-element  boundaries.  The  error  is  then  estimated  by  solving  a  local 
(i.e.,  element- by-element)  problem  using  either  the  principle  of  minimum  potential  en¬ 
ergy,  or  its  dual,  the  principle  of  minimum  complementary  energy,  on  a  space  which  is 
larger  than  the  one  used  for  obtaining  the  finite  element  solution.  The  space  on  which 
the  minimum  is  sought  can  be  enlarged  by  p-  or  fi-extension.  The  difficulty  with  achiev¬ 
ing  local  equilibrium  is  that  the  jump  discontinuities  must  be  carefully  distributed  such 
that  all  jumps  are  accounted  for  and  each  element  is  in  equilibrium  under  the  action 
of  the  distributed  load  corresponding  to  the  residuums  and  the  assigned  tractions  cor¬ 
responding  to  the  jump  discontinuities  in  stresses  along  the  perimeter. 

One  important  advantage  of  the  method  developed  in  this  project  is  that  equilibration 
is  not  required.  Recognizing  that  in  the  complementary  energy  principle  the  displace¬ 
ments  are  the  natural  boundary  conditions,  which  can  be  easily  applied,  the  local  prob¬ 
lem  is  formulated  using  the  displacements  computed  from  the  finite  element  solution 
and  a  stress  space  which  satisfies  both  the  equilibrium  and  compatibility  conditions. 
The  complementary  potential  energy  is  then  minimized  on  this  space.  The  method 
can  be  applied  on  individual  elements,  a  group  of  elements,  or  part  of  an  element. 
The  error  is  measured  by  the  work  of  the  jump  in  tractions  at  element  boundaries 
computed  from  the  solution  of  the  dual  problem,  using  the  displacement  computed 
from  the  finite  element  solution  (the  primal)  and  the  strain  energy  of  the  difference  in 
stresses  computed  from  the  primal  and  the  dual  on  the  element  domain.  The  sum  of 
the  error  measures  assigned  to  the  element  edges  and  to  the  elements  is  used  as  global 
error  estimator.  The  ratio  of  this  global  error  measure  and  the  total  potential  energy 
is  found  to  be  close  to  the  square  of  the  relative  error  in  energy  norm  for  a  wide  range 
of  problems.  Details  and  examples  are  available  in  [2]. 

2.  Model-adaptive  methods  for  laminated  plates  and  shells.  An  investigation  of 
model  adaptive  procedures  for  laminated  plates  and  shells  was  undertaken.  The  main 
points  are  outlined  in  a  highly  simplified  setting.  Details  are  available  in  [1]. 

Consider  an  elastic  strip  of  rectangular  cross-section  and  assume  that  the  x-y  plane 
is  a  principal  plane  of  the  cross  section  and  the  beams  are  subjected  to  bending  and 
shear  in  the  x-y  plane  only  and  the  axial  force  is  zero.  The  material  is  assumed  to  be 
elastic  and  isotropic.  Considering  this  as  a  two-dimensional  problem  of  elasticity,  we 
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can  write  the  displacement  vector  components  ux  and  uy  in  the  form: 


(3) 

(4) 


Ux  =  Uxll(x)4>i(y)  +  ux\3(x)My)  +  Ux\5(x)<t>b{y)  +  •••  +  U*|(2m-l)(z)02n»-l(y) 
Uy  =  Uy\o(x)lpo{y)  +  Uy]2(x)lp2(y)  +  Uy\i{x)ll)z{y )  H - h  UyK2n-2)(x)i’2n-2(y) 


where  <f>i  (resp.  pi)  are  antisymmetric  (resp.  symmetric)  director  functions  chosen  so 
that  an  optimal  rate  of  convergence  is  realized  in  the  sense  of  eq.  (7).  The  functions 
ux\ k(x),  uy\k(x )  are  called  field  functions. 

Writing  the  displacement  components  in  this  form  allows  one  to  consider  a  family  of 
hierarchic  models  (semi-discretizations)  for  beams  and  laminated  elastic  strips  char¬ 
acterized  by  m  and  n.  The  highest  member  of  the  hierarchy  is  the  problem  of  two- 
dimensional  elasticity  which  corresponds  to  m  =  n  =  oo.  Denoting  the  exact  solution 
of  the  hierarchic  model  characterized  by  m  and  n  by  an<4  the  exact  solution  of 

the  problem  of  two-dimensional  elasticity  by  the  hierarchic  models  satisfy  the 

following  three  criteria: 


(a)  Approximability: 


lim  Wu^x'1 


Uex  We 


=  0. 


where  1 1  ■  |  |jsr  is  the  energy  norm, 
(b)  Asymptotic  consistency: 


(5) 


lim 
h->  o 


\UEX 


—  U 


(m,n)  i 
EX  I 


,7(2D)  I 


U 


EX 


=  0. 


(6) 


(c)  Optimality  of  convergence: 


i?7(2D)  _  I 
I  uex  uex  I 


u 


(2  D), 
EX  I 


Ch?^rn'n'1  as  h i  —*  0,  m,  n  — >  oo 


(7) 


with  convergence  rates  7 (m  +  1,  n)  >  7 (m,  n);  7 (m,  n  +  1)  >  7(771,  n). 

The  formulation  and  analysis  of  hierarchic  models  for  plates  and  shells  was  completed. 
This  is  documented  in  [1]. 


3.  Asymptotic  expansions.  The  extraction  of  the  coefficients  of  the  eigenpairs  (the 
generalized  flux  intensity  factors  or  GFIFs  and  the  generalized  stress  intensity  fac¬ 
tors  or  GSIFs)  from  the  finite  element  solution,  at  multi-material  interfaces,  as  shown 
schematically  in  Fig.  1,  based  on  the  L2  and  energy  projection  method  on  a  circular 
sector,  was  investigated  [15].  The  problem  of  coupled  thermo-elastic  problems,  ac¬ 
counting  for  the  dependence  of  material  properties  on  temperature  was  investigated 


page  6/19 


and  implemented  for  problems  of  heat  conduction.  To  solve  this  nonlinear  problem,  in 
the  first  step  a  linear  heat  conduction  problem  is  solved  using  an  average  value  for  the 
coefficients  of  thermal  expansion.  In  subsequent  steps  the  thermal  conductivity  tensor 
is  determined  from  the  temperature  of  the  preceding  iteration  step  at  each  integra¬ 
tion  point  and,  using  these  values,  the  stiffness  matrix  is  re-computed.  The  stopping 
criterion  is  based  on  the  difference  of  two  consecutive  solution  vectors  [6],  [15]. 


Figure  1:  Typical  multi-material  interface. 


4.  Nonlinear  models:  Average  stress/strain  criteria.  Real  materials  typically  have 
a  periodic  internal  structure  characterized  by  the  grain  size,  molecular  length,  fiber 
diameters  and  the  spacing  of  the  fibers.  The  theory  of  elasticity  is  not  applicable  on 
length  scales  smaller  than  a  few  periods,  called  the  characteristic  length.  In  composite 
materials  it  is  often  necessary  to  account  for  material  nonlinearities,  however,  especially 
at  adhesively  bonded  joints.  For  these  reasons  the  use  of  the  average  stress  as  the  basis 
for  the  development  of  failure  criteria  is  far  more  realistic  than  the  determination  of 
the  first  term  of  the  asymptotic  expansion.  Nevertheless,  it  is  very  useful  to  have 
the  asymptotic  expansion  computed,  since  the  average  stress  on  a  small  area  can  be 
computed  very  efficiently  from  the  asymptotic  expansion. 

In  the  linear  case  it  can  be  seen  that  criteria  based  on  stress  or  strain  averages  are 
equivalent  to  the  criteria  based  on  stress  intensity  factors.  Specifically:  the  average 
stress  a(n,  A An)  is  defined  as: 


a(m,  A An) 


1 

ai; 


aijiriiUj  dA 


where  A An  is  a  small  area  the  orientation  of  which  is  characterized  by  the  unit  vector 
n,-.  The  direction  of  the  stress  acting  on  A An  in  given  by  the  unit  vector  m*,  as  shown 
in  Fig.  2. 
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Figure  2:  Notation. 

In  linear  elastic  fracture  mechanics  (LEFM)  this  is  equivalent  to  computing  the  stress 
intensity  factor.  For  example,  considering  a  cracked  plate  of  thickness  t  and  making 
the  usual  assumption  that  the  crack  is  aligned  with  the  aq-axis  and  the  crack  tip  is 
located  in  the  origin,  we  have: 

Ki  o  /  .  e  .  3 e\  _ .  1/2. 

1722  =  ~3Wr  C°S  2  (  +  Sm  2  S“  T  J  +  °(r 

Letting  m,  =  n*  =  {0  1}  (i.e.,  6  =  0)  and  A An  =  tA£  we  have 

tA£  J0  \p2/nr  y/2~Ai 

It  is  seen  that  the  average  stress  or  a  fixed  small  area,  characterized  in  this  example 
by  A£,  is  proportional  to  the  stress  intensity  factor  Kj.  The  proportionality  constant 
depends  on  AT  In  a  perfectly  homogeneous  and  isotropic  material  A£  is  an  arbitrary 
small  distance  selected  such  that  the  0( Ai1^2)  terms  are  negligible  in  comparison  with 
the  first  term.  Having  a  capability  for  computing  average  stresses  for  linear  as  well 
as  nonlinear  models  makes  it  possible  to  investigate  a  variety  of  failure  criteria  of  the 
type 


C  MAX  TmaX 

where  a  (resp.  f)  is  the  normal  (resp.  shear)  stress  on  a  given  small  area  [12],  [11],  [4], 
[5], 

In  addition,  the  problem  of  mechanical  contact  was  investigated  using  adaptive  meshing 
and  space  enrichment  methods  [7]. 
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5.  Limits  of  elastic  stability.  The  increasing  emphasis  on  affordability  in  military 
systems  has  led  to  a  number  of  advances  in  airframe  design  and  production.  Unitized 
airframe  structural  components  are  replacing  sheet  metal  built-up  structures  in  order 
to  reduce  part  count  and  assembly  cycle  times  and  costs.  The  qualification  of  thick 
aluminum  plate  stock  enables  the  machining  of  integral  structure  with  a  minimum 
of  expensive  tooling.  The  development  of  high  speed  machining  (HSM)  techniques 
has  further  enabled  the  fabrication  of  thin  lightweight  structure  to  provide  improved 
performance  along  with  affordability.  A  unitized  structural  component  is  shown  in 
Fig.  3.  The  web  thickness  is  in  the  range  2.5  to  5  mm,  the  flange  thickness  is  between 
2.5  and  3  mm.  The  development  of  methods  for  the  design  and  certification  of  unitized 
structural  components  based  on  numerical  simulation  is  of  obvious  importance.  One 
of  the  failure  modes  is  loss  of  stability. 


Figure  3:  A  unitized  structural  component. 

In  many  cases  failure  of  structural  components  is  caused  by  some  form  of  instability, 
such  as  buckling,  crippling  and  flutter.  Reliable  mathematical  models  that  account  for 
instability  will  make  it  possible  to  obtain  improved  estimates  of  the  margin  of  safety  of 
structural  components.  The  mathematical  models  used  in  current  professional  practice 
are  based  on  dimensionally  reduced  descriptions  of  an  elastic  body.  Such  models  cannot 
represent  the  initial  stress  state  where  the  assumptions  incorporated  in  dimensionally 
reduced  models  do  not  hold,  and  hence  may  lead  to  engineering  decisions  that  cause 
unexpected  failures  or  avoidable  weight  penalties. 
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An  investigation  of  linearized  models  of  structural  stability  was  completed.  The  for¬ 
mulation  of  the  problem,  documented  in  reference  [3],  requires  the  computation  of  the 
characteristic  values  of  non-compact  operators.  While  the  numerical  results  indicate 
that  a  broad  range  of  important  practical  problems  can  be  solved  by  this  method, 
questions  of  a  theoretical  nature  remained.  In  particular,  it  was  not  known  whether 
spurious  eigenvalues  would  appear  when  the  initial  stresses  are  bounded,  and  hew  sin¬ 
gular  points  might  pollute  the  computed  eigenvalues.  Of  practical,  engineering  interest 
are  problems  where  the  domain,  or  a  large  part  of  the  domain,  is  thin,  as  in  the  example 
shown  in  Fig.  3.  This  class  of  problems  is  of  particular  interest  in  aerospace  applica¬ 
tions  where  weight  considerations  dictate  efficient  use  of  material.  In  a  collaborative 
work  with  Professor  Manil  Suri  of  the  University  of  Maryland,  Baltimore  Count)’,  these 
problems  were  substantially  resolved. 

Two  mathematical  models  of  the  general  theory  of  elastic  stability  exist  in  the  classical 
literature:  These  models  are  called  the  Biot-Prager  model  and  the  Trefftz  model.  Both 
were  implemented  in  fully  three-dimensional  setting  so  that  numerical  experiments 
could  performed.  For  thin  structures  the  results  closely  matched  the  classical  results. 
It  was  found  that  the  two  classical  models  yield  virtually  identical  results. 

Some  fundamental  questions  concerning  the  existence  of  a  solution,  the  properties  of 
the  spectrum,  and  their  relationship  to  loss  of  stability  had  not  been  investigated.  In 
this  investigation  two  working  hypotheses  were  advanced  [3]: 

(a)  The  spectrum  is  a  point  spectrum,  hence  it  is  meaningful  to  consider  the  lowest 
nonzero  eigenvalue  as  an  indicator  of  the  onset  of  instability; 

(b)  The  minimal  eigenvalue  of  the  finite  dimensional  problem  converges  to  its  infinite 
dimensional  counterpart  as  the  finite  element  space  is  enlarged  (i.e.,  the  degrees 
of  freedom  are  increased). 

In  a  parallel  investigation  the  both  working  hypothesis  were  proven  by  Professors 
Manil  Suri,  Ivo  Babuska  and  Monique  Dauge.  The  second  hypothesis  was  shown  to 
hold  subject  to  the  condition  that  the  stress  field  is  bounded  and  the  solution  domain 
is  thin  [14].  In  many  practical  problems  this  condition  cannot  be  satisfied,  because 
inclusion  of  fillets,  needed  for  eliminating  stress  singularities  in  the  elastic  solution 

6.  Viscoelastic  models.  Investigation  of  the  problem  of  viscoelastic  solids  was  under¬ 
taken.  The  goal  is  to  represent  the  time-dependent  behavior  of  materials  using  the  p- 
and  hp- versions  of  the  finite  element  method  and  establish  a  robust,  reliable  method  for 
concurrent  control  of  the  spatial  and  temporal  errors.  Since  the  constitutive  equations 
are  of  the  form 

P  (D)<t  =  Q  (D)e 
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where  P (D)  and  Q (D)  are  differential  operators  in  time  and  cr  (resp.  e)  is  the  stress 
(resp.  strain)  tensor,  the  mixed  method  must  be  used.  The  mixed  formulation  in  the 
context  of  the  p-version  has  been  implemented  for  planar  and  axisvmmetric  prolems 
in  an  experimental  program  and  model  problems  have  been  solved. 

An  important  result,  obtained  recently  by  C.  Schwab  and  his  associates1  is  that  for  dif¬ 
fusion  problems  exponential  rate  of  convergence  can  be  realized  when  the  discontinuous 
Galerkin  method  (DG)  is  employed  utilizing  the  hp-version  in  space  and  time. 

In  viscoelastic  problems  and  problems  of  plasticity  the  stresses  are  not  computable  from 
the  displacement  field.  For  this  reason  it  is  necessary  to  employ  mixed  formulations, 
that  is,  formulations  in  which  the  displacement  and  stress  fields  are  approximated 
separately. 

First  we  consider  the  problem  of  linear  isotropic  elasticity.  The  solution  domain  is 
denoted  by  G.  The  body  force  vector  is  denoted  by  /*  and  the  tractions  acting  on  the 
boundary  T  is  denoted  by  <?,.  The  displacement  field  is  denoted  by  u  =  uf,  ( i  =  1,2). 
The  strain  tensor  is  defined  by: 


=  \  («y  +  ««)■  (8) 

The  stress  tensor  is  denoted  by  a  =  cr^-  and  the  formulation  is  stated  as  follows:  Find 
(u,  <r)  e  Hl(n)  x  i?°(G)  such  that  for  all  (v,w)  e  i^J(G)  x 


J  J  crijeJjtdS  =  J J  fiVitdS  +  J  QiVitdL 
J J  efjWijtdS  =  J  J  VijAjkiWkitdS 


(9) 

(10) 


where  Hk(Cl)  represents  the  usual  Sobolev  space  with  k  generalized  derivatives  on 
G,  with  L2(G)  =  H°(Q).  The  norm  of  Hk(Q)  will  be  denoted  by  ||  •  \\Hk.  Aijki  is 
a  symmetric  positive  definite  tensor2  which  depends  on  the  Young’s  modulus  E  and 
Poisson’s  ratio  v  and  whether  planes  stress  or  plane  strain  conditions  are  assumed,  t 
is  the  thickness. 


In  the  numerical  solution  we  use  the  spaces  14  =  14(A,p)  C  Ff^G),  Sh  =  Sh (A,  q)  C 
H°(Q)  where  A  is  the  finite  element  mesh  and  p  (resp.  q)  is  the  polynomial  degree 

'Schwab,  C.  and  Schotzau,  D.,  “Time  Discretization  of  Parabolic  Problems  by  the  hp-Version  of  the  Dis¬ 
continuous  Galerkin  Finite  Element  Method”,  Research  Report  No.  99-04  Seminar  fur  Angewandte  Mathe- 
matik,  ETH,  Zurich,  1999. 

2Aijki<Jijcrki  >  0  for  any  ^  0  in  H°(Q). 
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assigned  to  the  elements  for  displacements  (resp.  stresses).  Introducing  the  notation: 

b(<r,v)=  J J  aijeJjtdS 

G  ij  k\W  kltdS 

G(v)  =f  J j  fiVitdS  +  J  giVitdL 

the  (u,cr)  formulation  can  be  written  in  the  following  form:  Find  (u,  cr)  G  14  x  5/,  G 
H 1  x  H°  such  that 

a(cr,  w)  —  6(w,  u)  =0  V  w  G  Sh  (11) 

6(or,v)=G(v)  V  v  G  14.  (12) 

In  order  to  establish  that  the  system  of  equations  corresponding  to  (11)  and  (12)  is 
solvable,  it  is  necessary  to  show  that  the  Babuska-Brezzi  condition,  also  known  as  the 
inf-sup  condition,  is  satisfied: 

sup  >  C|| w||//o  V  w  G  Nh  C  H°(Q)  (13) 

cr eM>  llcrll H° 

sup  —  >  Cllvllffi  V  v  G  Vh  C  H'ifl)  (14) 

cresh  \\cr\\H0 

where  the  space  Nh  C  Sh  is  defined  by: 

Nh  =f  {cr  €  Sh  such  that  b(cr,  v)  =  0  V  v  G  14}. 

Noting  that  a(cr,  w)  is  coercive  on  H°  x  H°,  inequality  (13)  is  satisfied.  Inequality  (14) 

is  clearly  satisfied  if  we  can  select  a y  proportional  to  ejj.  For  quadrilateral  elements, 

product  space,  and  triangular  elements  this  is  possible  when  q  >  p—  1.  For  quadrilateral 
elements,  trunk  space3,  this  is  possible  if  q  >  p,  unless  the  trunk  space  is  so  constructed 
that  the  shape  functions  span  all  polynomials  of  degree  q  or  less,  but  not  more  than 
degree  q.  In  other  words,  the  xqy  and  xyq  terms  are  omitted,  in  which  case  q  >p—l. 
Omission  of  the  terms  xqy  and  xyq  is  possible  because  continuity  conditions  are  not 
enforced  on  Sh- 

Remark  1  Strictly  speaking,  the  rules  for  solvability  apply  to  linear  mapping  only.  In 
the  p  version,  where  generally  large  elements  are  used,  curved  elements  are  necessary 

3  Also  known  as  ‘serendipity’  space. 
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to  represent  curved  domains  accurately.  The  numerical  examples  presented  in  this 
paper  indicate  that  the  algebraic  equations  are  solvable  when  smooth  mapping  is  used. 
This  is  consistent  with  the  results  obtained  by  Chilton  and  Suri4  where  it  is  shown 
that  mixed  curved  quadrilateral  elements  retain  the  robustness  of  their  straight-sided 
counterparts,  provided  that  degenerate  elements  are  avoided. 


Summary 

The  problem  of  establishing  criteria  for  the  computation  of  margins  of  safety  for  laminated 
composites  and  adhesively  bonded  joints  requires  (a)  the  construction  of  a  mathematical 
model  that  accounts  for  the  essential  physical  attributes;  (b)  reliable  numerical  solution  of 
the  mathematical  model,  and  (c)  experimental  determination  of  the  allowable  values.  At 
present  the  principles  underlying  the  construction  of  mathematical  models  are  well  under¬ 
stood.  Efficient  methods  for  computing  the  data  of  interest  to  within  a  prescribed  tolerance 
are  now  available  commercially  and  research  is  being  performed  at  various  academic  insti¬ 
tutions  aimed  at  improving  the  efficiency  of  numerical  solution  methods.  A  great  deal  of 
work  remains  to  be  done  in  the  development  of  experimental  data  for  the  allowable  values, 
however.  The  experimental  work  must  be  coupled  with  reliable  numerical  solution  methods, 
the  foundations  for  which  have  been  laid  in  this  research  project. 
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Interactions/Transitions 

Industrial  application  has  been  made  possible  by  Engineering  Software  Research  and 
Development,  Inc.  (ESRD),  located  in  St.  Louis,  Missouri.  ESRD,  founded  in  1989  by 
Barna  Szabo,  Ivo  Babuska  and  Kent  Myers  develops,  documents  and  maintains  an  advanced 
professional  quality  finite  element  analysis  code  called  StressCheck  which  is  the  principal 
channel  by  which  the  results  of  AFOSR-sponsored  research  performed  at  the  Center  for 

4 Chilton  L,  Suri  M.  On  the  construction  of  stable  curvilinear  p  version  elements  for  mixed  formulations 
of  elasticity  and  Stokes’  flow.  Numer.  Math.  2000;  86:29-48. 
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Computational  Mechanics  of  Washington  University  have  been  transitioned  into  professional 
practice. 

The  current  users  of  StressCheck  include  The  Boeing  Aircraft  Company  (several  divi¬ 
sions);  Northrop  Grumman  Corporation;  Lockheed  Martin  Corporation;  Cessna  Aircraft 
Co.;  New  Piper  Aircraft  Co.,  Raytheon,  NASA  Johnson  Space  Center,  Israel  Aircraft  Indus¬ 
tries,  and  others. 

Specific  examples  of  recent  applications  by  persons  not  connected  to  this  project  are: 

1.  Report  by  C.  L.  Brooks  et  al.,  entitled  “Corrosion  is  a  Structural  and  Economic  Prob¬ 
lem:  Transforming  Metrics  to  a  Life  Prediction  Method”.  In  this  report  it  is  shown 
how  StressCheck  is  used  for  the  prediction  of  life  of  corroded  aircraft  components.  The 
work  described  in  the  report  was  performed  under  the  sponsorship  of  the  U.S.  Air 
Force. 

2.  Paper  by  Arnold  Nathan  of  the  Israel  Aircraft  Industries:.  “Composite  Repair  of  Aging 
Metallic  Structure:  p- Version  3D  Finite  Element  Approach”,  presented  at  the  2nd 
Joint  NASA/FAA/DoD  Conference  on  Aging  Aircraft,  31  August-3  September,  1998 
Williamsburg,  VA.  This  paper  describes  a  typical  application  of  StressCheck  to  the 
analysis  composite  material  repair  patch  bonded  to  a  metallic  structure. 

3.  Paper  S.  P.  Engelstad  and  R.  L.  Actis  [13]  The  desire  to  reduce  costs  of  composite 
aircraft  structure  has  recently  increased  the  importance  of  bonded  joint  analysis  ca¬ 
pabilities  in  the  aeronautics  industry.  A  new  initiative  has  been  formed,  called  the 
Composites  Affordability  Initiative  (CAI),  with  the  goal  to  decrease  the  recurring  ac¬ 
quisition  costs  of  composite  airframe  structure  by  50%,'  thereby  making  composites 
more  affordable  for  the  next  generation  of  fighter  aircraft.  Participants  in  the  initia¬ 
tive  include  the  Air  Force,  Navy,  Lockheed  Martin,  Boeing,  and  Northrop  Grumman. 
New  validated  analysis  tools  for  composite  bonded  joints  have  been  viewed  as  the  key 
enabler  for  the  application  of  innovative  manufacturing  and  bonding  technologies  for 
airframe  primary  structure.  Thus,  improved  analysis  tools  for  design  of  bonded  com¬ 
posite  joints  have  become  the  focus  of  a  CAI  analysis  tools  team.  This  paper  reviews 
the  development  history  and  describes  new  capabilities  of  a  p- version  Finite  Element- 
based  tool  developed  by  a  joint  effort  of  the  CAI  team  and  a  software  provider. 

The  capability  of  performing  rapid  parametric  sizing  studies  was  one  of  the  central  goals 
of  the  team.  It  was  desired  to  develop  parameterized  solutions  for  various  bonded  joint 
configurations,  and  validate  failure  predictions  with  CAI  test  data.  These  solutions 
would  then  be  used  to  perform  rapid,  preliminary  design  sizing  of  composite  bonded 
joints.  Bonded  joint  analysis  tools  previously  existed  in  the  aeronautics  industry  (such 
as  the  industry  standard  A4EI  code),  but  are  limited  to  simple  joint  geometry  (typ¬ 
ically  single,  double,  and  step-lap  joints)  due  to  the  one  dimensional  nature  of  the 
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solution.  Much  more  complex  joints  are  being  utilized  to  achieve  affordability  goals, 
thus  solution  techniques  require  at  a  minimum  two-dimensional  continuum  elasticity 
analysis.  Many  finite  element  studies  utilizing  plane  strain  analysis  appear  in  the  lit¬ 
erature,  all  utilizing  h-element  finite  element  models.  Many  utilized  nonlinear  models 
of  the  adhesive,  and  others  included  geometric  nonlinearity.  Different  methods  for 
handling  the  so-called  spew  fillet  geometry,  which  is  formed  at  the  adhesive  free-  edge, 
have  also  been  considered.  Uncertainty  of  the  shape  of  this  region  has  been  studied 
also.  Since  singularities  exist  in  the  elasticity  solution  at  free  edges  of  composite  layer 
interfaces,  and  also  at  reentrant  corners,  consideration  of  the  convergence  issues  and 
error  control  in  these  regions  is  important.  It  is  the  authors’  opinion  that  more  at¬ 
tention  to  error  control  in  these  regions  is  required,  since  failure  tends  to  initiate  in 
the  vicinity  of  these  singularities.  Other  important  issues  include  the  need  for  a  mod¬ 
ern,  user-friendly  interface,  abilities  for  the  user  to  build-in  failure  criteria,  advanced 
stress /strain  extraction  capabilities,  and  insensitivity  to  very  large  aspect  ratios  caused 
by  the  modeling  of  thin  composite  layers. 

The  CAI  analysis  tools  team  completed  an  extensive  study  of  the  state-of-the-art  in  the 
area  of  composite  stress  analysis  tools,  as  would  be  applied  to  failure  analysis  of  com¬ 
posite  bonded  joints.  The  software  StressCheck,  developed  by  Engineering  Software 
Research  and  Development,  Inc.  (ESRD),  was  identified  as  a  potential  provider,  and 
following  an  evaluation  phase,  was  selected.  This  software  offers  a  p-  and  hp-element 
formulation,  which  satisfies  the  need  for  error  control,  aspect  ratio  insensitivity,  and 
also  includes  geometric  and  material  nonlinearity.  StressCheck  already  had  an  elec¬ 
tronic  handbook  interface,  which  contained  many  of  the  required  parameterization 
features.  Once  StressCheck  became  an  approved  CAI  analysis  tool  for  further  devel¬ 
opment,  CAI  and  ESRD  worked  closely  to  define  StressCheck  modifications  required 
for  the  bonded  joint  applications.  ESRD  has  implemented  CAI  requested  software  im¬ 
provements  through  both  industry  contract  funding  and  AFOSR  sponsorship5.  This 
paper  presents  a  review  of  the  modifications  implemented  into  StressCheck,  and  a 
model  problem  consisting  of  a  double-lap  joint  analysis.  The  differences  .between  the 
new  and  conventional  modeling  procedures,  including  the  handling  of  singularities,  and 
error  control,  are  identified. 

This  paper  also  illustrates  an  example  of  a  very  fruitful  effort  of  the  CAI  team  in  iden¬ 
tifying  and  filling  an  analysis  tool  requirement.  An  important  conclusion  is  the  need 
for  industry  users  to  work  closely  with  software  providers  to  precisely  tailor  software 
to  the  requirements  of  the  problem  and  the  focus  user  group.  This  experience  high¬ 
lights  the  key  link  being  forged  between  aerospace  industry  researchers  and  software 
providers  to  develop  next  generation  analysis  tools. 

5Project  No.  FQ8671-9501469  STTR/TS. 
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The  Principal  Investigator  presented  the  results  obtained  under  this  project  to  defense  con¬ 
tractors  and  an  Air  Force  Laboratory: 

1.  The  Boeing  Company,  St.  Louis,  MO  (contact  person:  Mark  Holly  314-232-1405) 
August  18,  1998. 

2.  Materials  and  Manufacturing  Directorate,  AFRL  WPAFB  July  21,  1998  (contact  per¬ 
sons:  Dr.  Richard  B.  Hall  937-255-9097,  Dr.  Nicholas  Pagano  937-255-1138) 

3.  Lockheed  Martin,  Marietta,  GA.  Met  with  the  Composites  Affordability  Initiative 
(CAI)  Group  CAI  is  a  government/industry  consortium.  Members  represent  each  of 
the  major  US  aerospace  companies  and  the  Air  Force  and  Navy.  The  principal  inves¬ 
tigator  made  a  presentation  to  a  committee  of  CAI  charged  with  the  responsibility  of 
developing  procedures  for  the  reliable  and  efficient  numerical  simulation  of  structural 
components  fabricated  of  composite  materials.  Following  the  presentation  an  in-depth 
study  of  the  numerical  procedures  developed  under  AFOSR  sponsorship  at  the  Cen¬ 
ter  for  Computational  Mechanics  of  Washington  University  was  undertaken  by  CAI. 
(Contact  person:  Dr.  Steven  P.  Engelstad  770-494-9714)  March  3,  1998.  Frequent- 
follow-ups. 

4.  Presentation  to  The  Boeing  Company  St.  Louis:  “Engineering  decisions  based  on  com¬ 
puted  information:  Time,  quality  and  cost”  April  7,  2000  (contact  person:  Eric  Meyer 
314-232-4097) 

5.  Boeing  Aircraft  Company,  Seattle,  WA  “New  Directions,  Methods  and  Tools  for  Stress 
Analysis”  17  Oct.  2000 

6.  Boeing  Aerospace,  Long  Beach,  CA  “Quality  Assurance  in  Engineering  Computations. 
Applications  to  structural  design  and  analysis”,  October  16,  2000  (contact  person: 
Patrick  Goggin  563-593-7367). 

The  Principal  Investigator  delivered  lectures/presentations  at  the  following  professional  meet¬ 
ings: 

•  Szabo,  B.  and  Actis,  R.  L.,  “Hierarchic  Models  in  the  Numerical  Simulation  of  Mechan¬ 
ical  Systems”,  Keynote  Lecture,  4th  World  Congress  on  Computational  Mechanics, 
Buenos  Aires,  Argentina,  29  June  -  2  July  1998. 

•  Szabo,  B.  A.  and  Actis,  R.  L,  “Analysis  of  bonded  and  fastened  repairs  by  the  p- 
version  of  the  finite  element  method”,  The  Fourth  Joint  DoD/FAA/NASA  Conference 
on  Aging  Aircraft,  St.  Louis,  MO  May  15-18,  2000. 
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•  Seminars  at  three  German  Universities:  Munich  Stuttgart  and  Hannover  at  the  invi¬ 
tation  of  Professor  Erwin  Stein.  Seminar  titles;  “Numerical  Simulation  of  Mechanical 
Systems;  Problems  of  Reliability  and  Complexity”  and  “Model- Adaptive  Processes  in 
the  Numerical  Simulation  of  Mechanical  Systems”,  March  7-12,  2000. 

•  Invited  lecture:  “Quality  Assurance  in  the  Numerical  Simulation  of  Mechanical  Sys¬ 
tems”,  Mathematische  Analyse  von  FEM  fur  Probleme  in  der  Mechanik,  Mathematis- 
ches  Forschungsinstitut,  Oberwolfach,  Germany,  February  1999. 

•  Invited  paper:  “Mathematical  Models  of  Fastened  Structural  Connections”  SAE  Gen¬ 
eral,  Corporate  and  Regional  Aviation  Meeting,  Wichita,  KS,  20-22  April  1999 

•  Invited  minisymposium  presentation:  “Linear  models  for  estimating  the  limits  of  elastic 
stability”  1999  SIAM  Annual  Meeting,  Atlanta,  GA  May  1999. 

•  Invited  plenary  lecture:  “Standardization  of  Design  and  Analysis  Procedures”  Imple¬ 
mentation  Road  Map  Conference  1999,  Dearborn,  Michigan,  October  1999. 

•  Keynote  address:  “Hierarchic  modeling.  Concepts  and  implementation”,  Advanced 
Design  and  Analysis  Conference,  Tokyo,  Japan  December  16,  2000. 

•  Invited  paper:  “Quality  Assurance  in  the  Numerical  Simulation  of  Mechanical  Sys¬ 
tems”,  The  5th  International  Conference  on  Computational  Structures  Technology, 
Leuven,  Belgium,  September  2000. 

•  Contributed  paper:  “Analysis  of  Bonded  and  Fastened  Joints  by  the  p- Version  of  the 
Finite  Element  Method”  The  Fourth  Joint  DoD/FAA/NASA  Conference  on  Aging 
Aircraft,  St.  Louis  15-18  May,  2000 

•  Invited  presentation:  “New  FEA-Based  Technology  for  Integral  Structures”,  Aircraft 
Structural  Integrity  Program  Conference,  Integral  Structure  Technical  Exchange,  San 
Antonio,  Texas,  December  4,  2000 


New  discoveries,  inventions  and  patent  disclosures 

The  findings  of  this  research  are  essential  for  quantitative  assessment  of  the  strength  of 
homogeneous  and  heterogeneous  solid  bodies.  Failure  theories  require  experimental  correla¬ 
tion  and  confirmation  before  they  could  be  properly  called  discoveries,  however. 

Honors  /  Awards 

The  principal  investigator  received  an  honorary  doctorate  (doctor  honoris  causa,  1998). 
He  is  an  external  member  of  the  Hungarian  Academy  of  Sciences  and  Fellow  of  the  US 
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Association  for  Computational  Mechanics.  The  principal  investigator  serves  on  the  editorial 
boards  of  Computers  and  Structures,  and  Engineering  with  Computers,  and  on  the  scientific 
advisory  boards  of  the  First  M.I.T.  Conference  on  Computational  Fluid  and  Solid  Mechanics 
(June  2001). 
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